/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2013-2025 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "CMULES.H"
#include "fvcSurfaceIntegrate.H"
#include "localEulerDdtScheme.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

template<class RdeltaTType, class RhoType, class SpType>
void Foam::MULES::correct
(
    const RdeltaTType& rDeltaT,
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phiCorr,
    const SpType& Sp
)
{
    Info<< "MULES: Correcting " << psi.name() << endl;

    scalarField psiIf(psi.size(), 0);
    fvc::surfaceIntegrate(psiIf, phiCorr);

    psi.primitiveFieldRef() =
    (
        (rho.primitiveField()*rDeltaT - Sp.primitiveField())
       *psi.primitiveField()
      - psiIf
    )/(rho.primitiveField()*rDeltaT - Sp.primitiveField());

    psi.correctBoundaryConditions();
}


template<class RhoType>
void Foam::MULES::correct
(
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phiCorr
)
{
    correct(rho, psi, phiCorr, zeroField());
}


template<class RhoType, class SpType>
void Foam::MULES::correct
(
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phiCorr,
    const SpType& Sp
)
{
    const fvMesh& mesh = psi.mesh();

    if (fv::localEulerDdt::enabled(mesh))
    {
        const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
        correct(rDeltaT, rho, psi, phiCorr, Sp);
    }
    else
    {
        const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
        correct(rDeltaT, rho, psi, phiCorr, Sp);
    }
}


template<class RhoType, class PsiMaxType, class PsiMinType>
void Foam::MULES::correct
(
    const control& controls,
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phiBD,
    surfaceScalarField& phiCorr,
    const PsiMaxType& psiMax,
    const PsiMinType& psiMin
)
{
    correct
    (
        controls,
        rho,
        psi,
        phiBD,
        phiCorr,
        zeroField(),
        psiMax,
        psiMin
    );
}


template
<
    class RhoType,
    class SpType,
    class PsiMaxType,
    class PsiMinType
>
void Foam::MULES::correct
(
    const control& controls,
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phiBD,
    surfaceScalarField& phiCorr,
    const SpType& Sp,
    const PsiMaxType& psiMax,
    const PsiMinType& psiMin
)
{
    const fvMesh& mesh = psi.mesh();

    if (fv::localEulerDdt::enabled(mesh))
    {
        const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);

        limitCorr
        (
            controls,
            rDeltaT,
            rho,
            psi,
            phiBD,
            phiCorr,
            Sp,
            psiMax,
            psiMin
        );

        correct(rDeltaT, rho, psi, phiCorr, Sp);
    }
    else
    {
        const scalar rDeltaT = 1.0/mesh.time().deltaTValue();

        limitCorr
        (
            controls,
            rDeltaT,
            rho,
            psi,
            phiBD,
            phiCorr,
            Sp,
            psiMax,
            psiMin
        );

        correct(rDeltaT, rho, psi, phiCorr, Sp);
    }
}


template
<
    class RhoType,
    class SpType,
    class PsiMaxType,
    class PsiMinType
>
void Foam::MULES::correct
(
    const control& controls,
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phiBD,
    UPtrList<surfaceScalarField>& phiCorrs,
    const SpType& Sp,
    const PsiMaxType& psiMax,
    const PsiMinType& psiMin
)
{
    const fvMesh& mesh = psi.mesh();

    surfaceScalarField phiCorr(phiCorrs[0]);
    for (label i = 1; i<phiCorrs.size(); i++)
    {
        phiCorr += phiCorrs[i];
    }

    if (fv::localEulerDdt::enabled(mesh))
    {
        const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);

        limitCorr
        (
            controls,
            rDeltaT,
            rho,
            psi,
            phiBD,
            phiCorrs,
            phiCorr,
            Sp,
            psiMax,
            psiMin
        );

        correct(rDeltaT, rho, psi, phiCorr, Sp);
    }
    else
    {
        const scalar rDeltaT = 1.0/mesh.time().deltaTValue();

        limitCorr
        (
            controls,
            rDeltaT,
            rho,
            psi,
            phiBD,
            phiCorrs,
            phiCorr,
            Sp,
            psiMax,
            psiMin
        );

        correct(rDeltaT, rho, psi, phiCorr, Sp);
    }
}


template
<
    class RdeltaTType,
    class RhoType,
    class SpType,
    class PsiMaxType,
    class PsiMinType
>
void Foam::MULES::limitCorr
(
    const control& controls,
    const RdeltaTType& rDeltaT,
    const RhoType& rho,
    const volScalarField& psi,
    const surfaceScalarField& phiBD,
    surfaceScalarField& phiCorr,
    const SpType& Sp,
    const PsiMaxType& psiMax,
    const PsiMinType& psiMin
)
{
    const fvMesh& mesh = psi.mesh();

    surfaceScalarField lambda
    (
        IOobject
        (
            "lambda",
            mesh.time().name(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        mesh,
        dimensionedScalar(dimless, 1)
    );

    // Correction equation source
    scalarField SuCorr
    (
        mesh.Vsc()().primitiveField()
       *(rho.primitiveField()*rDeltaT - Sp.primitiveField())
       *psi.primitiveField()

    );

    limiter
    (
        controls,
        lambda,
        rDeltaT,
        rho,
        psi,
        SuCorr,
        phiBD,
        phiCorr,
        Sp,
        psiMax,
        psiMin
    );

    phiCorr *= lambda;
}


template
<
    class RdeltaTType,
    class RhoType,
    class SpType,
    class PsiMaxType,
    class PsiMinType
>
void Foam::MULES::limitCorr
(
    const control& controls,
    const RdeltaTType& rDeltaT,
    const RhoType& rho,
    const volScalarField& psi,
    const surfaceScalarField& phiBD,
    UPtrList<surfaceScalarField>& phiCorrs,
    surfaceScalarField& phiCorr,
    const SpType& Sp,
    const PsiMaxType& psiMax,
    const PsiMinType& psiMin
)
{
    const fvMesh& mesh = psi.mesh();

    surfaceScalarField lambda
    (
        IOobject
        (
            "lambda",
            mesh.time().name(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        mesh,
        dimensionedScalar(dimless, 1)
    );

    // Correction equation source
    scalarField SuCorr
    (
        mesh.Vsc()().primitiveField()
       *(rho.primitiveField()*rDeltaT - Sp.primitiveField())
       *psi.primitiveField()

    );

    limiter
    (
        controls,
        lambda,
        rDeltaT,
        rho,
        psi,
        SuCorr,
        phiBD,
        phiCorr,
        Sp,
        psiMax,
        psiMin
    );

    forAll(phiCorrs, i)
    {
        phiCorrs[i] *= lambda;
    }
    phiCorr *= lambda;
}


template
<
    class RhoType,
    class SpType,
    class PsiMaxType,
    class PsiMinType
>
void Foam::MULES::limitCorr
(
    const control& controls,
    const RhoType& rho,
    const volScalarField& psi,
    const surfaceScalarField& phiBD,
    surfaceScalarField& phiCorr,
    const SpType& Sp,
    const PsiMaxType& psiMax,
    const PsiMinType& psiMin
)
{
    const fvMesh& mesh = psi.mesh();

    if (fv::localEulerDdt::enabled(mesh))
    {
        const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);

        limitCorr
        (
            controls,
            rDeltaT,
            rho,
            psi,
            phiBD,
            phiCorr,
            Sp,
            psiMax,
            psiMin
        );
    }
    else
    {
        const scalar rDeltaT = 1.0/mesh.time().deltaTValue();

        limitCorr
        (
            controls,
            rDeltaT,
            rho,
            psi,
            phiBD,
            phiCorr,
            Sp,
            psiMax,
            psiMin
        );
    }
}


// ************************************************************************* //
